架空植物の第i番目の個体、個体ごとに、個体サイズ\(x_{i}\)と種子数\(y_{i}\)の依存関係を見る
個体差を想定しない場合には、
glm(y~x, family = poisson, data = ...)
これによって切片と傾きの最尤推定値が得られる
サイズによる種子数のバラツキは平均\(\lambda_{i}\)のポアソン分布\(p(y_{i}\mid \lambda_{i})\)に従うとする。
この例題では個体差をそうていしていないので、ランダム効果の項はなし。
ベイズモデルの事後分布は尤度×事前分布に比例するので、以下の様に表せる。
\[p(\beta_{1},\beta_{2}\mid \textbf{Y})\propto p(\textbf{Y}\mid \beta_{1},\beta_{2})p(\beta_{1})p(\beta_{2})\]
事前に確率分布がわからない場合には、線形予測子のパラメーター\(\beta\)が\([-\infty,\infty]\)の範囲で好きな値をとって良いという無情報事前分布を設定する
無情報っぽい事前分布として扱うのは、ばらつきがかなり大きくなる様に(偏りが出ない様に)
のいずれかを使う(この本では後者を使う)
winbugs使えないしstanを使う
インストール
インストールするのめっちゃめんどくせぇ。。。
インストール補足
catalinaからのインストールにはR4.0.0が最短かも
参考リンク
kubo9.stanファイルにパラメータなどベイズ統計モデルを記述
その後繰り返し回数と、burnin(最初の方の試行のうち、何回までを使わないか)を決める
なお、繰り返し回数については、複数のサンプル列(MCMCサンプル)を比較し、その大小によって収束しているかどうかを判断する
data {
int<lower=0> N;
real X[N];
real MeanX;
int<lower=0> y[N];
}
parameters {
real beta1;
real beta2;
}
model {
for (i in 1:N){
y[i] ~ poisson(exp(beta1 + beta2 * (X[i] - MeanX)));
}
beta1 ~ normal(0, 100); //無情報事前分布
beta2 ~ normal(0, 100); //無情報事前分布
}
library("rstan")
## Loading required package: StanHeaders
## Loading required package: ggplot2
## rstan (Version 2.19.3, GitRev: 2e1f913d3ca3)
## For execution on a local, multicore CPU with excess RAM we recommend calling
## options(mc.cores = parallel::detectCores()).
## To avoid recompilation of unchanged Stan programs, we recommend calling
## rstan_options(auto_write = TRUE)
load("/Users/hamaokigaizaburo/Documents/kuboR/intro-statistical-modeling-master/data/d.RData")
d.dat <- list(N=dim(d)[1], X=d$x, y=d$y, MeanX=mean(d$x)) #Stanに渡すデータ
d.fit <- stan(file='kubo9.stan', data=d.dat, iter=1600, warmup=100, thin=3, chains=3)
## Warning in readLines(file, warn = TRUE): '/Users/hamaokigaizaburo/Documents/
## kuboR/kubo9_files/kubo9.stan' で不完全な最終行が見つかりました
## Trying to compile a simple C file
d.fit
## Inference for Stan model: kubo9.
## 3 chains, each with iter=1600; warmup=100; thin=3;
## post-warmup draws per chain=500, total post-warmup draws=1500.
##
## mean se_mean sd 2.5% 25% 50% 75% 97.5% n_eff Rhat
## beta1 1.98 0.00 0.08 1.80 1.92 1.98 2.04 2.14 1456 1
## beta2 0.08 0.00 0.07 -0.05 0.04 0.08 0.13 0.23 1535 1
## lp__ 143.97 0.03 1.05 141.29 143.58 144.32 144.70 144.95 1326 1
##
## Samples were drawn using NUTS(diag_e) at Mon May 25 16:37:16 2020.
## For each parameter, n_eff is a crude measure of effective sample size,
## and Rhat is the potential scale reduction factor on split chains (at
## convergence, Rhat=1).
これにより、各パラメータの事後分布にまつわる情報を出すことができる
ex)beta1の値の範囲は95%の事後確率で、1.8 ~ 2.12 になる
また、収束指数はRhat列、有効なサンプルサイズはn_eff列に記述されている
このとき、隣り合うMCMCステップのサンプル間の相関が高い時に、この有効なサンプルサイズは小さくなる
plot(d.fit)
## ci_level: 0.8 (80% intervals)
## outer_level: 0.95 (95% intervals)
traceplot(d.fit)
library(ggplot2)
library(reshape2)
library(grid)
d.ext <- extract(d.fit, permuted=F)
b1.1 <- d.ext[1:500,'chain:1','beta1']
b1.2 <- d.ext[1:500,'chain:2','beta1']
b1.3 <- d.ext[1:500,'chain:3','beta1']
b1 <- data.frame(chain1=b1.1, chain2=b1.2, chain3=b1.3)
b1.melt <- melt(b1, id=c(), variable.name="chain") #ggplotで扱いやすいようにb1を再形成
b1.melt <- data.frame(b1.melt, iter=1:500)
#サンプリング過程
p <- ggplot(b1.melt, aes(x=iter, y=value, group=chain, color=chain))
p <- p + geom_line(size=0.1)
p <- p + labs(title="beta1のサンプリング過程", x="Iterations", y="")
#事後分布
g <- ggplot(b1.melt, aes(x=value))
g <- g + geom_density() + theme_bw()
g <- g + labs(title="beta1の事後分布", x="", y="")
#サンプリング過程と事後分布を並べて表示
grid.newpage()
pushViewport(viewport(layout=grid.layout(1, 2)))
print(p, vp=viewport(layout.pos.row=1, layout.pos.col=1))
print(g, vp=viewport(layout.pos.row=1, layout.pos.col=2))
MCMCサンプリングのためのアルゴリズム
この新しい値の確率分布とは、多変量確率分布からひとつの変量をのぞいて、他の変量すべてを定数とする一変量確率分布。これを全条件つき分布という